Load required packages

library(tidyverse)
library(phyloseq)
library(speedyseq)
library(ggrepel)
library(ampvis2)
library(plotly)
library(microbiome)
library(here)
options(getClass.msg=FALSE) # https://github.com/epurdom/clusterExperiment/issues/66
#this fixes an error message that pops up because the class 'Annotated' is defined in two different packages

Source required functions

rm(list = ls())

'%!in%' <- function(x,y)!('%in%'(x,y))

source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_taxa_tests.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_normalisation.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_alpha.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_beta.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")

Import phyloseq object

ps = "data/raw/metabarcoding/ps_silva_dada2_human_chicken_meta.RDS"

ps %>% 
  here::here() %>%
  readRDS() %>%
  phyloseq_get_strains_fast() %>% 
  filter_taxa(function(x) sum(x > 0) > 0, TRUE) %>% 
  subset_samples(Enrichment == "NotEnriched") %>% 
  subset_samples(Model == "Human")  %>% 
    subset_samples(Treatment %!in% c("positive", "negative", "DONOR", "CCUG59168", "HV292.1", "STRAIN")) -> physeq
## Joining, by = "ASV"
physeq %>% 
  sample_data() %>% 
  data.frame() %>% 
  DT::datatable()

Rarefaction:

physeq %>% 
  rarefy_even_depth(sample.size = 4576,
                    rngseed = 123) -> physeq_rare
## `set.seed(123)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(123); .Random.seed` for the full vector
## ...
## 38 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## IR-13-S103IR-46-S126IR-51-S157IR-74-S188IR-82-S143
## ...
## 738OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
physeq_rare
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 372 taxa and 123 samples ]:
## sample_data() Sample Data:        [ 123 samples by 59 sample variables ]:
## tax_table()   Taxonomy Table:     [ 372 taxa by 8 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 372 tips and 371 internal nodes ]:
## refseq()      DNAStringSet     :      [ 372 reference sequences ]
## taxa are rows
physeq_rare %>% 
  phyloseq_compute_bdiv(phylo = FALSE, norm = "pc") -> bdiv_list
## Loading required package: ape
## Loading required package: GUniFrac
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
## 
##     diversity
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
physeq_rare  %>%
  phyloseq_plot_bdiv(bdiv_list,
                     m = "CoDa", # with this potion I am generating CLR transformed OTU table
                     seed = 123) -> coda

bdiv_list$coda <- coda$physeq_clr %>% phyloseq::distance(method = "euclidean")  # Aitchinson distance wich claim to be superior to the other ones for compositional data.
coda$PCA$layers[[1]] = NULL

coda$PCA + geom_point(size=2,
                      aes(color =Treatment, 
                          fill = Treatment,
                          shape = NULL,
                          alpha = Day_of_Treatment)) + 
  theme_light() +
  geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.15, "inches")),
            size = 0.05, linetype = "dashed", inherit.aes = TRUE, aes(group=Reactor_Treatment, color = Treatment), show.legend = FALSE) +
  scale_alpha_continuous(range=c( 0.9, 0.3)) + scale_color_viridis_d(na.value = "black") + 
  scale_fill_viridis_d(na.value = "black") + 
  # scale_shape_manual(values = c(8, 21, 22, 23, 24, 16, 15, 18, 17)) + 
  theme_classic() +
  labs(col=NULL, fill = NULL, shape = NULL) + guides(shape=FALSE) -> pca_atchi

pca_atchi

plotly::ggplotly(pca_atchi)
# pca_atchi %>% 
#   export::graph2ppt(append = TRUE,
#                     file = file.path(here::here("data/processed/figures_NRP72")))

This is the fonction of interest

phyloseq_add_taxa_vector(dist =  bdiv_list$coda,
                         phyloseq = physeq_rare,
                         figure_ord = pca_atchi,
                         pval_cutoff = 0.05,
                         top_r = 16, # top features tp display(strain, or depending on what you specified based on correlation r value)
                         taxrank_glom = "Strain",
                         tax_rank_plot = "Strain"
) -> out

Add the vectors to the ordination plot

out$plot

Get the raw table of taxa and correlation with the ordination:

out$envfit %>%
  DT::datatable()

You could then filter that table and use the ASV id to generate a heatmap of only the significant features:

out$envfit %>%
  filter(pval <= 0.05) %>% 
  pull(id) -> ASV_envit
physeq_rare %>% 
prune_taxa(x = ., taxa = ASV_envit) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Model", "Fermentation", "Reactor_Treatment"),
                          tax_aggregate = "Genus",
                          tax_add = NULL,
                          ntax  = 50) -> p

p$data %>% 
  mutate(Abundance = na_if(Abundance, 0)) -> p$data 


p + facet_grid( ~Model + Fermentation + Reactor_Treatment , scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 0.1, 0.2, 0.50, 0.75, 0.100), 
                       labels = c(0,  0.01, 0.1, 0.2, 0.50, 0.75, 0.100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2

# p2 %>% 
  # export::graph2ppt(append = TRUE,
  #                   file = file.path(here::here("data/processed/figures_NRP72")))
vegan::adonis(bdiv_list$coda  ~ get_variable(physeq_rare, "Reactor_Treatment") * 
                get_variable(physeq_rare, "Day_of_Treatment"))$aov.tab %>% 
  data.frame() %>%
  rownames_to_column('Term') %>%
  mutate_if(is.numeric, round, 4) %>%
  DT::datatable()
lapply(
  bdiv_list,
  FUN = phyloseq_adonis_strata_perm,
  physeq = physeq_rare,
  formula = paste0(c("Reactor_Treatment", "Day_of_Treatment"), collapse=" * "),
  nrep = 999,
  strata = "none"
) %>%
  bind_rows(.id = "Distance") %>%
  mutate_if(is.numeric, round, 3) %>%
  DT::datatable()
phyloseq_plot_bdiv(physeq_rare,
                   bdiv_list,
                   m = "NMDS",
                   axis1 = 1,
                   axis2 = 2) -> plots
# plot_list %>%
#   phyloseq_ordinations_expl_var()
physeq_rare %>%
  phyloseq_distance_boxplot(bdiv_list$coda ,
                            d = "Reactor_Treatment") -> dist_box


dist_box_plot <- dist_box$plot

dist_box_plot

# 
# dist_box_plot %>% 
#   export::graph2ppt(append = TRUE,
#                     file = file.path(here::here("data/processed/figures_NRP72")))
plots %>%
  phyloseq_plot_ordinations_facet(color_group =  "Day_of_Treatment",
                                  shape_group = "Treatment")  -> plot_list

plot_list

lapply(
  bdiv_list,
  FUN = physeq_pairwise_permanovas,
  physeq = physeq_rare,
  compare_header = "Reactor_Treatment",
  n_perm = 999,
  strat = FALSE
) %>%
  bind_rows(.id = "Distance") %>%
  mutate_if(is.numeric, round, 3) %>%
  # filter(! terms %in% (c("Residuals", "Total"))) %>%
  DT::datatable()
correct_plot_x_continious <- function(original_plot = tmp_plot$data){
  
  
  ggplot(data=original_plot,mapping=aes(x=as.numeric(as.character(varGroup2)),y=Distance,color=Label)) +
    # geom_boxplot(data=tmp_plot$data,mapping=aes(x=as.factor(varGroup2),y=Distance,color=Label),outlier.size = 0.5) +
    geom_point(position=position_jitterdodge(jitter.width = 0.1,seed=123),aes(group=Label),size=2, alpha=0.4) +
    theme_bw() + 
    geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")),
                           size = 0.4, linetype = "dashed", inherit.aes = TRUE, aes(group=Label, color = Label),
                           position=position_jitterdodge(dodge.width=0.9)) +
    theme(legend.position="bottom")  + geom_smooth(show.legend = FALSE, level = 0.95, alpha=0.005, size = 0.001) + 
    # ylim(c(0,1)) +
    xlab("Day Treatment") + guides(col = guide_legend(ncol = 3)) + scale_x_continuous(breaks=seq(0,90,10)) + scale_fill_viridis_d() + scale_color_viridis_d() -> plot
  
  return(plot)
}
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda", "bjaccard", "wjaccard"),
                                     bdiv_list = bdiv_list,
                                     physeq = physeq_rare,
                                     timepoint = "fixed",
                                     group_var = "Reactor_Treatment",
                                     time_var = "Day_of_Treatment",
                                     fixed_time = -2) -> t
tmp_plot <- t$bjaccard + facet_null() 

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected 

# tmp_plot_corrected %>% 
#   export::graph2ppt(append = TRUE,
#                     file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$wjaccard + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected

# 
# tmp_plot_corrected %>% 
#   export::graph2ppt(append = TRUE,
#                     file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$coda + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected %>% 
  plotly::ggplotly()
# tmp_plot_corrected %>% 
#   export::graph2ppt(append = TRUE,
#                     file = file.path(here::here("data/processed/figures_NRP72")))
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda","bjaccard", "wjaccard"),
                                     bdiv_list = bdiv_list,
                                     physeq = physeq_rare,
                                     timepoint = "previous",
                                     group_var = "Reactor_Treatment",
                                     time_var = "Day_of_Treatment") -> t


tmp_plot <- t$bjaccard + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot <- t$wjaccard + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot <- t$coda + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) -> tmp_plot_corrected

tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda","bjaccard", "wjaccard"),
                                     bdiv_list = bdiv_list,
                                     physeq = physeq_rare,
                                     timepoint = "between.ref.group",
                                     group_var = "Reactor_Treatment",
                                     time_var = "Day_of_Treatment",
                                     group_to_compare = "CR_UNTREATED") -> t
tmp_plot<- t$bjaccard + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected

tmp_plot_corrected

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$wjaccard + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected

tmp_plot_corrected

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$coda + facet_null()

tmp_plot$data %>% 
  correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected

tmp_plot_corrected

tmp_plot_corrected %>% 
  export::graph2ppt(append = TRUE,
                    file = file.path(here::here("data/processed/figures_NRP72")))
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gdtools_0.2.2        usedist_0.4.0.9000   GUniFrac_1.1        
##  [4] matrixStats_0.58.0   vegan_2.5-7          lattice_0.20-41     
##  [7] permute_0.9-5        ape_5.4-1            reshape2_1.4.4      
## [10] scales_1.1.1         here_1.0.1           microbiome_1.10.0   
## [13] plotly_4.9.2.1       ampvis2_2.6.5        ggrepel_0.8.2       
## [16] speedyseq_0.5.0      phyloseq_1.34.0      forcats_0.5.0       
## [19] stringr_1.4.0        dplyr_1.0.4          purrr_0.3.4         
## [22] readr_1.4.0          tidyr_1.1.2          tibble_3.0.6        
## [25] ggplot2_3.3.3        tidyverse_1.3.0.9000
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1            uuid_0.1-4              backports_1.2.1        
##   [4] systemfonts_0.3.2       plyr_1.8.6              igraph_1.2.6           
##   [7] lazyeval_0.2.2          splines_4.0.2           crosstalk_1.1.0.1      
##  [10] digest_0.6.27           foreach_1.5.1           htmltools_0.5.1.1      
##  [13] magrittr_2.0.1          cluster_2.1.0           doParallel_1.0.16      
##  [16] openxlsx_4.2.3          Biostrings_2.56.0       modelr_0.1.8           
##  [19] officer_0.3.14          prettyunits_1.1.1       colorspace_2.0-0       
##  [22] blob_1.2.1              rvest_0.3.6             haven_2.3.1            
##  [25] xfun_0.21               crayon_1.4.1            jsonlite_1.7.2         
##  [28] survival_3.2-3          iterators_1.0.13        glue_1.4.2             
##  [31] rvg_0.2.5               gtable_0.3.0            zlibbioc_1.34.0        
##  [34] XVector_0.28.0          webshot_0.5.2           Rhdf5lib_1.10.1        
##  [37] BiocGenerics_0.34.0     DBI_1.1.1               miniUI_0.1.1.1         
##  [40] Rcpp_1.0.6              xtable_1.8-4            viridisLite_0.3.0      
##  [43] progress_1.2.2          stats4_4.0.2            DT_0.15                
##  [46] truncnorm_1.0-8         htmlwidgets_1.5.3       httr_1.4.2             
##  [49] RColorBrewer_1.1-2      ellipsis_0.3.1          pkgconfig_2.0.3        
##  [52] NADA_1.6-1.1            farver_2.0.3            dbplyr_1.4.4           
##  [55] manipulateWidget_0.10.1 later_1.1.0.1           tidyselect_1.1.0       
##  [58] labeling_0.4.2          rlang_0.4.10            munsell_0.5.0          
##  [61] cellranger_1.1.0        tools_4.0.2             cli_2.3.0              
##  [64] generics_0.1.0          devEMF_4.0-2            ade4_1.7-16            
##  [67] export_0.3.0            broom_0.7.2             fastmap_1.1.0          
##  [70] evaluate_0.14           biomformat_1.7.0        yaml_2.2.1             
##  [73] knitr_1.31              fs_1.5.0                zip_2.1.1              
##  [76] rgl_0.100.54            nlme_3.1-149            mime_0.10              
##  [79] xml2_1.3.2              compiler_4.0.2          rstudioapi_0.13        
##  [82] zCompositions_1.3.4     reprex_0.3.0            stringi_1.5.3          
##  [85] highr_0.8               stargazer_5.2.2         Matrix_1.2-18          
##  [88] multtest_2.44.0         vctrs_0.3.6             pillar_1.4.7           
##  [91] lifecycle_1.0.0         data.table_1.13.6       flextable_0.5.11       
##  [94] httpuv_1.5.4            patchwork_1.1.0         R6_2.5.0               
##  [97] promises_1.1.1          network_1.16.0          IRanges_2.22.2         
## [100] codetools_0.2-16        ggnet_0.1.0             MASS_7.3-52            
## [103] assertthat_0.2.1        rhdf5_2.32.4            rprojroot_2.0.2        
## [106] withr_2.4.1             S4Vectors_0.26.1        mgcv_1.8-32            
## [109] parallel_4.0.2          hms_1.0.0               grid_4.0.2             
## [112] rmarkdown_2.4           Cairo_1.5-12.2          Rtsne_0.15             
## [115] shiny_1.5.0             Biobase_2.50.0          lubridate_1.7.9        
## [118] base64enc_0.1-3